library(janitor)library(readxl) # to read in an .xlsx filelibrary(car) library(rstanarm)library(easystats)library(tidyverse)theme_set(theme_bw())knitr::opts_chunk$set(comment =NA)source("c08/data/Love-431.R") # for the lovedist() function
Importing Data / Creating Factors
ohio20 <-read_xlsx("c08/data/ohio_2020.xlsx") |>mutate(behavior = Hmisc::cut2(rk_behavior, g =4),clin_care = Hmisc::cut2(rk_clin_care, g =3)) |>mutate(behavior =fct_recode(behavior,"Best"="[ 1,23)", "High"="[23,45)","Low"="[45,67)", "Worst"="[67,88]")) |>mutate(clin_care =fct_recode(clin_care,"Strong"="[ 1,31)", "Middle"="[31,60)","Weak"="[60,88]")) |>select(FIPS, state, county, outcomes, behavior, clin_care, everything())
Can we identify pairwise comparisons of means that show meaningful differences?
Can we build a confidence interval procedure that protects against Type I error expansion due to multiple comparisons?
Two Methods for Multiple Comparisons
There are two (main) methods we’ll study to estimate the difference between specific pairs of means with uncertainty intervals, while dealing with the problem of multiple comparisons.
If we do 6 tests, we could reduce the necessary \(\alpha\) to 0.01 / 6 = 0.00167 and that maintains an error rate no higher than \(\alpha = 0.01\) across the 6 tests.
Better Approach: Holm-Bonferroni
Suppose you have \(m\) comparisons, with p-values sorted from low to high as \(p_1\), \(p_2\), …, \(p_m\).
Is \(p_1 < \alpha/m\)? If so, reject \(H_1\) and continue, otherwise STOP.
Is \(p_2 < \alpha/(m-1)\)? If so, reject \(H_2\) and continue, else STOP.
and so on…
Holm-Bonferroni Approach
This is uniformly more powerful than Bonferroni, while preserving the overall false positive rate at \(\alpha\). It’s also the default choice of estimate_contrasts().
estimate_contrasts(model_one, ci =0.99, p_adjust ="holm")
# A tibble: 6 × 10
Level1 Level2 Difference CI_low CI_high SE df t p contr
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
1 Best High 0.711 0.204 1.22 0.156 84 4.55 3.52e- 5 Best - Hi…
2 Best Low 1.14 0.635 1.65 0.156 84 7.32 5.53e-10 Best - Low
3 Best Worst 2.01 1.50 2.51 0.156 84 12.8 9.54e-21 Best - Wo…
4 High Low 0.431 -0.0759 0.939 0.156 84 2.76 7.04e- 3 High - Low
5 High Worst 1.29 0.787 1.80 0.156 84 8.29 7.87e-12 High - Wo…
6 Low Worst 0.863 0.356 1.37 0.156 84 5.53 1.07e- 6 Low - Wor…
Plot the Holm results
ggplot(con_holm_tib, aes(x = contr, y = Difference)) +geom_point(size =3, col ="purple") +geom_errorbar(aes(ymin = CI_low, ymax = CI_high)) +geom_hline(yintercept =0, col ="red", lty ="dashed") +labs(title ="Holm 99% Intervals for Health Outcomes",x ="Contrast", y ="Difference in Health Outcomes Score")
Plot the Holm results
Tukey’s HSD Intervals
Tukey’s HSD approach is a better choice for pre-planned comparisons with a balanced (or nearly balanced) design.
estimate_contrasts(model_one, ci =0.99, p_adjust ="tukey")
con_ols <-as_tibble(estimate_contrasts(model_one, ci =0.99, p_adjust ="Holm")) |>mutate(contr =str_c(Level1, " - ", Level2))con_bay <-as_tibble(estimate_contrasts(model_1b, ci =0.99, p_adjust ="Holm")) |>mutate(contr =str_c(Level1, " - ", Level2))p1 <-ggplot(con_ols, aes(y = contr, x = Difference)) +geom_point(size =3, color ="firebrick") +geom_errorbar(aes(xmin = CI_low, xmax = CI_high)) +geom_vline(xintercept =0, col ="red", lty ="dashed") +labs(title ="OLS fit",subtitle ="Holm 99% HSD Intervals for Health Outcomes",y ="Contrast", x ="Difference in Health Outcomes Score")p2 <-ggplot(con_bay, aes(y = contr, x = Difference)) +geom_point(size =3, color ="navy") +geom_errorbar(aes(xmin = CI_low, xmax = CI_high)) +geom_vline(xintercept =0, col ="red", lty ="dashed") +labs(title ="Bayesian fit",subtitle ="Holm 99% HSD Intervals for Health Outcomes",y ="Contrast", x ="Difference in Health Outcomes Score")p1 + p2
Comparison of Holm 99% CIs
K Samples: Comparing Means
What is the outcome under study?
What are the (in this case, \(K \geq 2\)) treatment/exposure groups?
Were the data in fact collected using independent samples?
Are the data random samples from the population(s) of interest? Or is there at least a reasonable argument for generalizing from the samples to the population(s)?
What is the confidence level we require?
Are we doing one-sided or two-sided testing? (usually 2-sided)
What does the distribution of each individual sample tell us about which inferential procedure to use?
Are there meaningful differences between population means?
Can we identify pairwise comparisons of means that show meaningful differences using an appropriate procedure that protects against Type I error expansion due to multiple comparisons?
